function [ X1, X2, P1, P2, BI1, BI2 ] = PspatialMapPerBlock( h , NIt , NJt , NJe , NIe , PT)
%SPATIALMAP Summary of this function goes here
%   Detailed explanation goes here

x1 = zeros (NIe*NJe,1);
y1 = zeros (NIe*NJe,1);
P1 = zeros (NIe*NJe,size(PT,2));

x2 = zeros ((NIt-NIe)*(NJt-NJe),1);
y2 = zeros ((NIt-NIe)*(NJt-NJe),1);
P2 = zeros ((NIt-NIe)*(NJt-NJe),size(PT,2));

xStart = h/2;
yStart = h/2;

index1 = 1;
index2 = 1;

for i = 1:NIt
    if(i<=NIe+1)
        x1(index1) = xStart + (i-1) * h;
        y1(index1) = 0;
        for t = 1:size(PT,2)
            P1(index1,t) = interp1([h,2*h, 3*h],[PT(Pmap(i,3),t),PT(Pmap(i,2),t), PT(Pmap(i,1),t)],[4*h],'cubic','extrap');
        end
        index1 = index1 + 1;
        x1(index1) = xStart + (i-1) * h;
        y1(index1) = yStart + (NJe-1/2) * h;
        for t = 1:size(PT,2)
            P1(index1,t) = interp1([3*h,2*h, h],[PT(Pmap(i,NJe-2),t),PT(Pmap(i,NJe-1),t), PT(Pmap(i,NJe),t)],[0],'cubic','extrap');
        end
        index1 = index1 + 1;
    else
        x2(index2) = xStart + (i-1) * h;
        y2(index2) = 0;
        for t = 1:size(PT,2)
            P2(index2,t) = interp1([h,2*h, 3*h],[PT(Pmap(i,3),t),PT(Pmap(i,2),t), PT(Pmap(i,1),t)],[4*h],'cubic','extrap');
        end
        index2 = index2 + 1;
        x2(index2) = xStart + (i-1) * h;
        y2(index2) = yStart + (NJt-1/2) * h;
        for t = 1:size(PT,2)
            P2(index2,t) = interp1([3*h,2*h, h],[PT(Pmap(i,NJt-2),t),PT(Pmap(i,NJt-1),t), PT(Pmap(i,NJt),t)],[0],'cubic','extrap');
        end
        index2 = index2 + 1;
    end
    if(i==NIe+1)
        x1(index1) = xStart + (i-1) * h;
        y1(index1) = 0;
        for t = 1:size(PT,2)
            P1(index1,t) = interp1([h,2*h, 3*h],[PT(Pmap(i,3),t),PT(Pmap(i,2),t), PT(Pmap(i,1),t)],[4*h],'cubic','extrap');
        end
        index1 = index1 + 1;
        x1(index1) = xStart + (i-1) * h;
        y1(index1) = yStart + (NJe-1/2) * h;
        for t = 1:size(PT,2)
            P1(index1,t) = interp1([3*h,2*h, h],[PT(Pmap(i,NJe-2),t),PT(Pmap(i,NJe-1),t), PT(Pmap(i,NJe),t)],[0],'cubic','extrap');
        end
        index1 = index1 + 1;
        x2(index2) = xStart + (i-1) * h;
        y2(index2) = 0;
        for t = 1:size(PT,2)
            P2(index2,t) = interp1([h,2*h, 3*h],[PT(Pmap(i,3),t),PT(Pmap(i,2),t), PT(Pmap(i,1),t)],[4*h],'cubic','extrap');
        end
        index2 = index2 + 1;
        x2(index2) = xStart + (i-1) * h;
        y2(index2) = yStart + (NJt-1/2) * h;
        for t = 1:size(PT,2)
            P2(index2,t) = interp1([3*h,2*h, h],[PT(Pmap(i,NJt-2),t),PT(Pmap(i,NJt-1),t), PT(Pmap(i,NJt),t)],[0],'cubic','extrap');
        end
        index2 = index2 + 1;
    end
end

BI1 = index1;
BI2 = index2;

for i=1:NIt
    for j=1:NJt
        
        if (i<=NIe && j>NJe) continue; end
        
        index = Pmap (i,j);
        
        if(i<=NIe)
            x1 (index1) = xStart + (i-1) * h; 
            y1 (index1) = yStart + (j-1) * h;
            P1(index1,:) = PT(index,:); 
            index1 = index1 + 1;
        else
            x2 (index2) = xStart + (i-1) * h; 
            y2 (index2) = yStart + (j-1) * h;
            P2(index2,:) = PT(index,:); 
            index2 = index2 + 1;
        end
        if(i==NIe+1 && j<=NJe)
            x1 (index1) = xStart + (i-1) * h; 
            y1 (index1) = yStart + (j-1) * h;
            P1(index1,:) = PT(index,:); 
            index1 = index1 + 1;            
        end
    end
end

X1 = [x1,y1];
X2 = [x2,y2];

end

